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Aggregation of stiff polyelectrolytes in solution 
and angle- and distance-dependent potential of 
mean force between two like-charged rods are 
studied in the presence of 3-valent salt using 
molecular dynamics simulations. In the bulk so- 
lution, formation of long-lived metastable struc- 
tures with similarities to the raft-like structures 
of actin filaments is observed within a range of 
salt concentration. The system finally goes to a 
state with lower free energy in which finite-sized 
bundles of parallel polyelectrolytes form. Pre- 
ferred angle and interaction type between two 
like-charged rods at different separations and 
salt concentrations are also studied, which shed 
some light on the formation of orientationally 
ordered structures. 



1 Introduction 

It is well known that highly charged anionic 
biological polyelectrolytes such as DNA and 
filamentous actin (F-actin) could attract each 
other in the presence of multivalent cations [1] . 



Attraction between similarly charged polyelec- 
trolytes, which is due to electrostatic correla- 
tions and can not be described in the frame- 
work of mean field theories such as Poisson- 
Boltzmann formalism, has been the subject of 
intense ongoing investigations in the past few 

decades |Ill2lEliai5lEllZlEllSllIDlinilI2l|T31. 

Aggregation of stiff polyelectrolytes in solution 
in the presence of multivalent salt arises from 
such like-charge attractions. 

Depending on their different roles in the cell 
function, actin filaments assemble into net- 
works (large crossing angle) or bundles (nearly 
parallel), respectively, helped by bundling or 
cross-linking proteins [14] . A similar pattern — 
namely, formation of both networks and dense 
bundles — has also been observed experimen- 
tally in solutions of F-actin rods with multiva- 
lent salt, depending on the salt concentration 
[T5] . In a particular range of salt concentra- 
tions, a multi-axial liquid crystalline phase of 
actin rods has been observed in which aggre- 
gates resembling stacks of raft-like structures 
form. These structures, which consist of a com- 
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bination of mutually parallel and perpendicular 
F-actin rods, change to bundles of parallel rods 
with increasing multivalent salt concentration 
[T5] . This is a manifestation of the rich com- 
plexity of high density aggregates of linear poly- 
electrolytes due to frustration in accommodat- 
ing strong electrostatic interactions and geomet- 
rical constraints simultaneously Similar effects 
arising from electrostatic frustration have been 
studied theoretically in the context of rodlike 
polyelectrolyte aggregates P2] , rodlike polyelec- 
trolyte brushes [IT] , and star polyelectrolytes 

Formation of novel raft-like aggregates of actin 
rods in the presence of multivalent salt has been 
studied both theoretically and using computer 
simulation methods. As a model system, ne- 
matic mixture of rigid rods and strong | cross- 
linkers has been studied theoretically using On- 
sager excluded volume theory, and the possibil- 
ity of formation of stacks of raft-like structures 
has been shown [19J. Distance-dependent inter- 
action between two parallel similarly charged 
rigid rods in the presence of explicit mono- 
and multivalent salt ions and angle-dependent 
interaction between them at close center-to- 
center separation have been studied by Lee et 
al. using molecular dynamics (MD) simulations 
[2D] . In their study, Lee et al. have shown 
that just above the threshold value of multiva- 
lent salt concentration needed for attraction be- 
tween parallel rods, the preferred angle at close 
separations is 90°, and they have concluded that 
this is in agreement with formation of raft-like 
structures in multivalent salt solution of F-actin 
rods. 

The phase behavior of a model system com- 
posed of similarly charged rigid rods and inter- 
rod linkers has also been studied using Debye- 
Hiickel theory of electrostatic interaction be- 
tween the rods and considering the role of the 
linkers in the effective inter-rod interaction po- 
tential [21]. In this study, it has been shown 
that aggregates with different orientational or- 
derings of the rods can be formed at different 
values of rod and linker concentrations. Re- 
cently, aggregation of stiff polyelectrolytes has 
been studied using MD simulations in the pres- 
ence of multivalent counterions [22J and in the 



multivalent salt solution [23]. In Ref. [23J it is 
shown that the dominant kinetic mode in aggre- 
gation process is the case of one end of one rod 
meeting others at right angle with no apprecia- 
ble energy barrier and the polyelectrolytes have 
a tendency to form finite-sized bundles. The ki- 
netics of aggregation and bundle formation of 
actin has also been recently probed experimen- 
tally [24]. Using two different fluorescently la- 
beled populations of F-actin, it was shown that 
the growth mode of actin bundles in the late 
stages is predominantly longitudinal, and that 
the bundles can freely exchange filaments with 
the solution and therefore are not in a frozen 
state [21]. It was also found that the energy 
barrier for the aggregation process is negligibly 
small (of the order of k^T) [21] in agreement 
with the theoretical result of Ref. ]23j . 

Considering the complexity of the problem of 
highly charged rod-like polyelectrolytes at high 
densities, it seems that more studies are needed 
to fully unravel the dynamic behavior of the sys- 
tem. Here we use a combination of simulation 
of the bulk solution of stiff polyelectrolytes and 
calculation of the effective interaction potential 
between a pair of polyelectrolytes in the pres- 
ence of multivalent salt to better understand the 
orientational ordering of the polyelectrolytes in 
the aggregation process. We use MD simula- 
tions to study stiff polyelectrolytes in the bulk 
solution in the presence of explicit counter-ions 
and ions of multivalent salt. We also examine 
the angle- and distance-dependent potential of 
mean force between two similarly charged rigid 
rods in the presence of multivalent salt to help 
understand the observations of the bulk solution 
simulations. We find that there could be three 
different regimes in the aggregation behavior of 
stiff polyelectrolytes in a multivalent salt solu- 
tion, depending the salt concentration. When 
the salt concentration is low, the system does 
not undergo aggregation because the attraction 
between polyelectrolytes is not strong enough 
to trigger it. In an intermediate range of salt 
concentrations, relatively long-lived structures 
with special orientational ordering of polyelec- 
trolytes form. In these scaffold-like structures, 
polyelectrolytes are predominantly joined either 
parallel or perpendicular to each other (see Fig. 
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[JJ), which is presumably due to the dominant 
kinetic mode of aggregation [23]. After suffi- 
ciently long time, however, mutually perpendic- 
ular rods in each structure slide against each 
other and join up in parallel, such that even- 
tually bundles of parallel polyelectrolytes form. 
These meta-stable structures have similarities 
to the raft-like structures of actin filaments ob- 
served experimentally [15]. In the high salt con- 
centration regime, the polyelectrolytes directly 
aggregate into bundles of parallel rods and do 
not exhibit right-angle configurations in their 
kinetic pathway of aggregation. 
We also calculate the angle- and distance- 
dependent potential of mean force between two 
similarly charged rigid rods, and find that in the 
low salt concentration regime there is no attrac- 
tion between rods. In the regime with interme- 
diate values of salt concentration, the interac- 
tion between parallel rods is attractive but in a 
range of center-to-center separation where the 
preferred angle between the rods is 90°. In this 
regime, the parallel configuration is preferred 
only at very close separations. In the high salt 
concentration regime, the parallel configuration 
is preferred at all separations and the interac- 
tion between parallel rods is attractive. 
We also check the stability (life-time) of raft- 
like structures by constructing such structures 
and monitoring their subsequent evolution as 
shown in Fig. [21 at different values of salt con- 
centration. We find that these structures are 
not stable in the low and high salt concentra- 
tion regimes: while at low salt concentrations 
mutual repulsion between rods destroys such 
structures, at high salt concentrations mutually 
perpendicular rods slide against each other in 
a very short time and the raft-like structures 
change to bundles of parallel rods. In the in- 
termediate salt concentration regime, raft-like 
structures are relatively long-lived. We have 
observed that the life time of such structures 
rapidly increases with the number of rods. 
The rest of the paper is organized as follows. In 
Sec. |2l the simulation method and the results 
of our studies on the aggregation of polyelec- 
trolytes are explained. Section [3] is devoted to 
the studies on the potential of mean force be- 
tween rods at various separations and angles, 



while Sec. @] describes the studies on the stabil- 
ity of raft-like structures. Finally, Sec. [5] con- 
cludes the paper. 

2 Aggregation of stiff poly- 
electrolytes in the bulk 
solution 

2.1 The model and the simula- 
tion method 

In our simulations, which are performed with 
MD simulation package ESPResSo v. 1.8 [2"5] . 
each polyelectrolyte rod is composed of N m = 
21 spherical charged monomers of charge — e 
(e is electronic charge) and diameter a. We 
have considered N p polyelectrolytes with their 
monomers being bonded to each other via a 
FENE [26] (finite extensible non linear elastic) 
potential, with a the separation between them 
being fixed at a = 1.1a. The bending rigid- 
ity of polyelectrolyte chains is modeled with 
bond angle potential = k^l — cos0) with 
kfj, = in which <fr is the angle between 

two successive bond vectors along the polyelec- 
trolyte chain. We use N c = N p x N m mono- 
valent counter-ions with charge +e to neutral- 
ize the polyelectrolyte charges. We also model 
3-valent salt as N s+ positive ions with charge 
+3e and N s _ = 3N S+ negative ions with charge 
— e. We follow Ref. [20] for the definition of 
salt concentration: c^.\ = N s _/N c (for exam- 
ple C3 : i = 0.5 means that iV s _ = 0.5 x N c ). 
The monomer number density in our simula- 
tions is 0.04 <t~ 3 , which is about two times the 
concentration we have used in Ref. [23J. In 
addition to long-range Coulomb interaction, we 
include short-range Lennard- Jones repulsion be- 
tween particles, which adds an energy scale e 
and a length scale a to the system. MD time 
step in our simulations is r = 0.01 tq in which 
t = J In ^- is the MD time scale and m is the 
mass of the particles. The temperature is fixed 
at k B T = 1.2e using a Langevin thermostat. We 
use the P3M method to apply periodic bound- 
ary conditions for long-range Coulomb interac- 
tion in the system. The strength of the elec- 
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trostatic interaction energy relative to the ther- 
mal energy can be quantified by definition of the 

2 

Bjerrum length /b = ^-f, where e is the dielec- 
tric constant of the solvent and in our simula- 
tions /b — 3.2a. In water at room temperature 
Zb = 7 A and the value of /b = 3.2a in our simu- 
lations means that a = 2.2A and the separation 
between charged monomers is a = 2.5A. By 
using the relation for the effective viscosity of 
the solvent i] ~ 2A^fmtj a 2 (from Ref. |27J), we 
can deduce an estimate for the microscopic time 
r ~ r/a 3 1 \2k^T) . Using the viscosity of water 
and room temperature, we find r = 1.3 ps. 
In the beginning of bulk solution simulations, 
we fix the polyelectrolytes in space and leave 
all the ions to fluctuate around them for 100000 
MD time steps, which is enough to equilibrate 
the ions for every given configuration of the 
polyelectrolytes. After equilibration, we remove 
the constraint on the polyelectrolytes and study 
the aggregation process in the system. 

2.2 Results 

Our bulk solution simulations show that with 
changing multivalent salt concentration, three 
different regimes can exist. At salt concentra- 
tions less than a lower value c^,\ ~ 0.2 no aggre- 
gation is observed in the system. In this regime, 
the salt concentration is lower than its mini- 
mum value needed for like-charge attraction be- 
tween polyelectrolytes. When a 3-valent salt ion 
links two monomers from two different polyelec- 
trolytes to each other, thermal fluctuations sep- 
arate them again and the system does not un- 
dergo an aggregation process. At intermediate 
values of the salt concentration 0.2 < c^i < 1.2, 
aggregation of polyelectrolytes can be observed 
and the dominant kinetic mode of aggregation 
is that of one end of one polyelectrolyte meeting 
others at a right angle, as previously reported 
|23j . During the aggregation process in this 
regime, long-lived structures of polyelectrolytes 
form, and in each structure polyelectrolytes are 
dominantly joined either parallel or perpendic- 
ular to each other (see Fig. H]f). 
In this regime, after sufficiently long time mu- 
tually perpendicular polyelectrolytes in each 
scaffold-like structure slide against each other 



and finally form bundles of parallel polyelec- 
trolytes. The kinetic mode of evolution from 
scaffold-like structures to bundles is consider- 
ably slower than the typical kinetic modes of 
the system. For example, in a system contain- 
ing nine polyelectrolytes, the time that takes for 
the scaffold-like structure (Fig. [1] f ) to evolve 
into the final configuration of a collapsed paral- 
lel bundle is greater than the time interval be- 
tween snapshots a and e of Fig. [TJ by a factor of 
at least six. This suggests that the formation of 
long-lived scaffold-like structures corresponds to 
meta-stable states of the system and the states 
with bundles of parallel polyelectrolytes have 
lower free energies. We have checked that the 
formation of the meta-stable structures is inde- 
pendent of the random number sequence used 
in the simulation, as well as the initial config- 
uration of the polyelectrolytes. In other words, 
in the intermediate salt concentration range the 
relaxation pathway to equilibrium appears to 
always involve intermediate meta-stable struc- 
tures. These structures have similarities to the 
raft-like structures of actin filaments observed 
experimentally [15] (see below for a discussion). 
In the regime with high values of salt concen- 
tration (c3 : i > 1.2), we observe a different ki- 
netic mode of aggregation in the system. In 
this regime, the polyelectrolytes directly aggre- 
gate into bundles of parallel polyelectrolytes and 
the salt concentration washes out the right angle 
configuration from the kinetic pathway of aggre- 
gation. The aggregation process in this regime 
is considerably faster than that of the regime 
with intermediate salt concentration. 
To shed some light on the behavior of the system 
and better understand these results, we study 
the potential of mean force between two poly- 
electrolyte rods. 

3 Potential of mean force 
between charged rods 

3.1 The model and the simula- 
tion method 

In our simulations of two similarly charged rigid 
rods we consider two rods, each composed of 
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N m = 21 similarly charged spherical monomers 
of diameter a and charge — e, N c = 2N m coun- 
terions of charge +e, and 3-valent salt of con- 
centration C3 : i. The monomer number density 
in these simulations is 10~ 4 a~ 3 , which corre- 
sponds to a simulation box of length L 75a 
that is large enough to avoid finite size effects. 
The rods are fixed on two parallel plates of sep- 
aration R and the angle between the rods is 9 
(see the schematic configuration in the inset of 
the middle part of Fig. [3]). For each set of 
values of R, 9, and c^-.i, the two rods are fixed 
during the simulation and only the counterions 
and the salt ions are free to move. After equi- 
libration of the system we calculate the aver- 
age force on each monomer and obtain the total 
force and the total torque around the center- 
to-center line (y-axis) for one of the rods. We 
integrate the force with respect to the center-to- 
center distance, R, to obtain the independence 
of the potential of mean force, and similarly in- 
tegrate the torque with respect to the angle 9 to 
obtain its ^-dependence (corresponding to rota- 
tion around y-axis). We set the zero of the po- 
tential of mean force at R = 7.5a for obtaining 
its independence, and at 9 = 90° in the calcu- 
lation of its ^-dependence. 

3.2 Results 

We first repeat the calculation of the potential 
of mean force between two parallel rods at dif- 
ferent values of salt concentration and center- 
to-center separation as in Ref. [20j. Our results 
show that for 3-valent salt concentrations higher 
than C3;i = 0.2, there is attraction between par- 
allel rods (see Fig. HJ), in agreement with the 
results of Ref. [20J. We now consider the fol- 
lowing question: if one allows the angle between 
the rods to vary for other values of R, will the 
rods prefer to be parallel? To answer this ques- 
tion, we calculate the angle-dependent potential 
of mean force at different values of R. We find 
that at the lowest value of R in our simulations, 
R = 2.5a, the preferred angle at salt concen- 
trations less than a threshold value c^-i — 0.2 is 
6 = 90°, and at higher salt concentrations, the 
parallel configuration is preferred (see R = 2.5a 
part of Fig. [3]). While this is the same result as 



reported in Ref. [20J , we find that upon increas- 
ing R (from R = 2.5a), the threshold value of 
salt concentration below which 9 = 90° is pre- 
ferred increases (see the parts of Fig. [3] corre- 
sponding to R = 3.5a and R = 4.5a ). These re- 
sults show that in a range of salt concentrations, 
if we consider two rod-like polyelectrolytes in 
the solution at a separation where their effec- 
tive interaction is appreciable, the interaction 
potential tends to reorient them to a perpen- 
dicular configuration. The preferred orientation 
is parallel only at very close separations. This 
could be the origin of the formation of scaffold- 
like structures in the aggregation of stiff poly- 
electrolytes. After the formation of such struc- 
tures, the polyelectrolytes are much closer to 
each other, and will have a chance to discover 
that the parallel configuration has a lower free 
energy, so that perpendicular rods eventually re- 
orient themselves and form bundles. 

To elaborate further on this issue, we calculate 
the distance-dependent potential of mean force 
between two perpendicular rods, which shows 
that in the presence of 3-valent salt there can 
be attractive interaction between them at close 
separations (see Fig. EJ). Moreover, Fig. [5] 
shows that a barrier exists in the potential of 
mean force, which decreases with increasing salt 
concentration. In Ref. [23J, it has been shown 
that the perpendicular configuration of two rods 
when their ends are close to each other has a 
lower free energy than the configuration shown 
in Fig. In fact, within the dominant kinetic 
mode of aggregation, when two polyelectrolytes 
are going to join each other with their ends and 
at right angle, they experience the lowest possi- 
ble energy barrier. Using the simulation results 
presented in Figs. El IU and El and similar stud- 
ies at other salt concentrations (not presented 
here), we have calculated the preferred crossing 
angle, 9 P , and the interaction type between two 
rods at different salt concentrations and center- 
to-center separations as shown in Fig. [HI We 
find that with increasing salt concentration, the 
9 = | region of the diagram becomes smaller, 
until it finally vanishes at c^ : \ ~ 1.0. Moreover, 
there is a region in the c$-i — R plane for which 
the preferred angle between the rods is 90° and 
the interaction is attractive. The existence of 
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such a regime could explain the formation of the 
scaffold-like structures during the aggregation 
of polyelectrolytes. We note that the informa- 
tion extracted from the potential of mean force 
between a pair of rods on the preferred angle 
between the rods and whether the interaction is 
attractive or repulsive does not provide a com- 
plete description of the behavior of the system, 
as the true phase behavior needs to be stud- 
ied from the global stability (convexity) of the 
free energy of the system (see e.g. Ref. [28] for 
such a study in the context of columnar DNA- 
aggregates). Such a study, however, will involve 
making assumptions on the structure of the ag- 
gregate as studying the full infinite dimensional 
configuration space of a system of N p rods is not 
computationally feasible. Instead, we prefer to 
use the calculation of the potential of the mean 
force (which is more tractable) to supplement 
the results obtained from MD simulations that 
probe the kinetics of the system while relaxing 
to equilibrium. 



We also calculate the ensemble average of inter- 
action energy (Lennard- Jones + Coulomb) and 
the potential of mean force (by integration of 
the torque with respect to the angle) as func- 
tions of the angle between rods with monova- 
lent counterions and no added salt. Potential of 
mean force AU is a partial free energy (subject 
to the constraint that the rods are kept fixed) 
and can be decomposed into an energetic and an 
entropic part, namely AU = AE — TAS, where 
the energetic contribution is the average inter- 
action energy AE. The potential of mean force 
AU and the average of interaction energy AE 
are shown in Fig. The comparison shows 
that although the right angle configuration is 
the preferred form judging from the potential 
of mean force, the parallel configuration has a 
lower interaction energy. This could be under- 
stood by noting that the counterions in the par- 
allel configuration are effectively confined in a 
smaller volume around the rods than in the per- 
pendicular configuration, which causes an en- 
tropic penalty. This suggests that the entropy 
of the counterions has a crucial role in determin- 
ing the equilibrium configuration of the system. 



4 Stability of the raft-like 
structures 

The formation of scaffold-like meta-stable struc- 
tures during the aggregation of rod-like poly- 
electrolytes has similarities to the raft-like struc- 
tures observed experimentally in solutions of 
F-actin rods with multivalent salt [15J. Since 
we observe that in an intermediate range of 
salt concentrations the system is likely to be 
trapped in meta-stable states on its way to equi- 
librium, the question naturally arises as to how 
this finding might have relevance to experiments 
on highly charged rod-like polyelectrolytes such 
as F-actin in the presence of multivalent salt. 
The formation of aggregates from single fila- 
ments distributed randomly in the solution as 
the initial configuration is controlled by a time 
scale that has two components, namely the 
time that takes the filaments to find each other 
through diffusion and the time that takes to 
overcome energetic barriers |23j . Since the diffu- 
sion component of the characteristic relaxation 
time could be substantial for dilute solutions, 
having long relaxation times does not necessar- 
ily mean that large energetic barriers impede 
the formation of equilibrium structures [23l [2^] . 
To eliminate the diffusion component of the ag- 
gregation time and focus specifically on the en- 
ergetic barriers in the dynamics of the system, 
we can study the evolution of already made 
complexes towards equilibrium. While strictly 
speaking this corresponds to a different kinetic 
pathway, one can imagine that a similar process 
of disentanglement of the interacting filaments 
will be followed in the aggregation kinetics, on 
the way to equilibrium. Therefore, measuring 
the relaxations times of already formed struc- 
tures could provide a direct quantitative mea- 
sure of the kinetics of the meta-stable state of 
the system. To this end, we study systemati- 
cally the stability and time evolution of raft-like 
aggregates under various conditions. 

We construct raft-like structures of stiff poly- 
electrolytes at different salt concentrations, and 
fix them at the beginning of the simulation (see 
Fig. [2] a) while the counterions and the salt ions 
are allowed to fluctuate. After the equilibration 
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of the free ions, we release the constraints on 
the polyelectrolyte rods and follow the evolu- 
tion of the system. We find that in the low salt 
concentration regime thermal fluctuations de- 
stroy these structures in a short time. In fact, 
in this regime, there is not enough salt in the 
system for like-charge attraction between poly- 
electrolytes to appear and the electrostatic re- 
pulsion between the rods make such structures 
unstable. In the high salt concentration regime 
also the raft-like structures have a fast evolution 
and considerably short life time. In this regime 
mutually perpendicular rods slide against each 
other in a short time and the structure changes 
to bundles of parallel rods. 
In the regime with intermediate salt concentra- 
tion, however, the system has a different evo- 
lution. In this regime we find that the raft- 
like structures are considerably long-lived and 
it takes a long time for the fluctuations to de- 
stroy them. The system finally finds the lower 
free energy state in which bundles of parallel 
polyelectrolytes are formed. In this regime, the 
stability of raft-like structures increases with in- 
creasing salt concentration. A time lapse view 
of the system with N p — 9 rods is shown in Fig. 
® 

To help quantify the stability of the raft-like 
structures, we study the dependence of the evo- 
lution time of the structures towards bundles 
of parallel polyelectrolytes, clS cl function of the 
size of these structures while keeping the vol- 
ume density of the filaments constant (by in- 
creasing the size of the simulation box with in- 
creasing number of filaments). The evolution 
time t e is defined as the time it takes for a raft- 
like structure containing N p polyelectrolytes to 
change to bundles of polyelectrolytes; for exam- 
ple the time interval between snapshots a and f 
of Fig. [2j The evolution time t e SIS db function 
of the number of polyelectrolytes in the raft- 
like structure N p is shown in Fig. [8^1. The 
data points in this figure are obtained by av- 
eraging over five different runs of the system. 
In this figure, N p = 9 corresponds to three mu- 
tually perpendicular layers and each layer con- 
tains three parallel rods (see Fig. [2]). Similarly, 
N p = 15, 20, 25, 30, 40 respectively correspond 
to 3, 4, 5, 6 and 8 layers and each layer contains 



five parallel rods. 

The characteristic evolution time t e can be 
thought of as the time it takes for the raft- 
like structure to disentangle the filaments that 
strongly interact with each other via frustrated 
electrostatic interactions, and as can be seen in 
Fig. [Hh. it increases with increasing the size 
of the raft-like structure N p . The increase is 
faster than linear, and can be approximated 
as a power law with an exponent of about 2.7 
(see the inset of Fig. [Hk for a log-log plot 
of t e versus N p ), although the numerical data 
are not over a sufficiently large domain so that 
the power law nature can be verified. The 
characteristic relaxation time of the system as 
a function of the number of filaments (for a 
fixed volume density of single filaments) can be 
used to estimate the N p dependence of the en- 
ergy barrier (or activation energy) E a for the 
transition between the raft-like structure and 
the equilibrium structure. Assuming a form of 
t e (N p ) = t e Ea( - Np ^ kBT , where to is a microscopic 
time scale that does not depend on the number 
of filaments, we can calculate relative changes 
in the activation energy as shown in Fig. [8b. 
While the plot in Fig. [Hb clearly shows that 
the energy barrier increases with the number of 
filaments, more studies are necessary to deter- 
mine whether this trend will continue or ulti- 
mately level off for systems with considerably 
larger number of filaments. 

Using the estimate for the microscopic time 
scale To = 1.3 ps, we find that for 9 rods the 
relaxation time is ~ 1 ns and for 40 rods it is 
~ 0.1 fxs. To make a crude estimate for the typ- 
ical relaxation times corresponding to realistic 
experimental cases [To], we can extrapolate the 
empirical scaling law to N p ~ 10 7 , which yields 
t e ~ 10 7 s. This estimate suggests that the time 
the system could be trapped in a meta-stable 
state could be considerably long compared to 
typical observations times, for a system like F- 
actin. Further studies are necessary to deter- 
mine whether the super-linear increase of the 
relaxation time as a function of the number of 
filaments persists for such large numbers present 
in the experiments on F-actin. 
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5 Concluding Remarks 

We have studied the aggregation of stiff poly- 
electrolytes in bulk solution in the presence of 
multivalent salt, and have observed spontaneous 
formation of orientationally ordered aggregates. 
Our results show that in a range of salt con- 
centrations, when two polyelectrolyte rods are 
at separations where their effective interaction 
is considerable, they could prefer to be per- 
pendicular while their interaction is attractive. 
In this range of salt concentrations, we have 
observed the formation of relatively long-lived 
scaffold-like structures of the polyelectrolytes in 
the bulk solution, for which the mode of evo- 
lution through to the bundles of parallel poly- 
electrolytes is considerably slower than the typ- 
ical kinetic modes of the system. We have also 
observed that in the above-mentioned range of 
the salt concentrations already formed raft-like 
structures similar to the ones observed in Ref. 
[T5] have longer life time than those in the high 
and low salt concentration regimes. 

Our findings suggest that when multivalent salt 
concentration is enough for like-charge attrac- 
tion between polyelectrolytes but is not so high 
that they join up in parallel configuration di- 
rectly, formation of raft-like structures could be 
possible as a meta-stable state. In our model 
system, which contains a small number of poly- 
electrolytes, the transition from the meta-stable 
state to the lower free energy state of bundles 
of parallel polyelectrolytes occur over a finite 
observable time. This evolution time, however, 
is observed to increase rapidly with the system 
size. This means that in a real experimental 
system with a large number polyelectrolytes in- 
volved in such a raft-like structure, it will be 
extremely difficult for the system to find the 
state with lower free energy in which the poly- 
electrolytes form dense bundles of parallel poly- 
electrolytes. In other words, our studies sug- 
gest that the experimentally observed raft-like 
structures might be meta-stable states with life 
times that are far larger than the experimen- 
tally available observation time. Further stud- 
ies on larger systems and over longer simulation 
times are needed to fully verify this proposition, 
and characterize the nature of the intermediate 



scaffold-like regime. 

It is a great pleasure to acknowledge stimulating 
discussions with Gerard Wong. 
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Fig. 1 Snapshots of the system contain- 
ing N p = 9 polyelectrolytes in so- 
lution with 3-valent salt concentra- 
tion c 3: i = 0.5 function of 
time. +3 salt ions are shown by 
golden (light) spheres, +1 counter- 
ions by red (dark), and -1 salt ions by 
gray spheres. Corresponding times 
of snapshots in units of r : (a) 10, 
(b) 25, (c) 40, (d) 55, (e) 70, (f) 85. 
Aggregates of polyelectrolytes with 
scaffold-like structures can be seen 
in the snapshots e and f. Finally, 
all of the polyelectrolytes form a sin- 
gle bundle (last configuration is not 
shown here). 

Fig. 2 Snapshots from evolution of a raft- 
like structure formed by N p — 9 
polyelectrolytes in the presence of 3- 
valent salt of concentration c^_\ = 
0.5. The positive 3-valent salt ions 
are shown by golden (light) spheres 
and the negative salt ions and coun- 
terions are not shown for clarity. 
Corresponding times of snapshots in 
units of r : (a) initial configuration 
before equilibration of the ions, (b) 
50, (c) 250, (d) 400, (e) 550, (f) 700. 

Fig. 3 Potential of mean force between two 
rods cLS cL function of 9 at different 
values of R and salt concentration 
C3 : i (see the schematic configuration 
of the rods in the inset of the mid- 
dle plot). The legends of the other 
two plots are the same as that of the 
rightmost plot. 

Fig. 4 Potential of mean force between two 
parallel rods (see the schematic con- 
figuration of the rods in the inset) as 
a function of center-to-center separa- 
tion R for different concentrations of 
3-valent salt. It can be seen that at 
salt concentrations c^.\ > 0.2 there is 
attraction between two parallel poly- 
electrolytes. 

Fig. 5 Potential of mean force between two 



perpendicular rods (schematic con- 
figuration of the rods is shown in the 
inset) cLS cL function of R at different 
concentrations of 3-valent salt. The 
size of the circles on the c^i = 
curve shows the size of error bars 
for all our data points. Note that 
the potential is not shown for R < 
1.6cr where the excluded volume in- 
teraction takes over and generates a 
strong repulsion. 

Fig. 6 Preferred angle 9 P and interaction 
type between two charged rigid rods 
which are fixed on two parallel plates 
of separation R, in different regions 
of C3;i — R plane (see the schematic 
configuration of the rods in the mid- 
dle plot of Fig. E]). 

Fig. 7 Ensemble average of the interaction 
energy (Lennard- Jones + Coulomb) 
per particle between two rods with 
their separation fixed at R = 2.5a as 
a function of 9 with no added salt. 
Inset: potential of mean force be- 
tween the same two rods. The dif- 
ferent trends of the two curves shows 
that the entropy of counterions has 
an important role in the orientation 
of the rods. 

Fig. 8 (a) Evolution time of raft-like struc- 
tures to bundles of parallel polyelec- 
trolytes as a function of N p , num- 
ber of polyelectrolytes forming these 
structures. Error bars are obtained 
from averaging over five runs of 
the simulation with different random 
number sequences. Inset: Log-log 
plot of t e versus N p and a line with 
the slope of 2.7 (dashed line), (b) 
Activation energy extracted from (a) 
as a function of N v . 
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Figure 1: Snapshots of the system containing 
N p = 9 polyelectrolytes in solution with 3-valent 
salt concentration c^-i = 0.5 function of 
time. +3 salt ions are shown by golden (light) 
spheres, +1 counter-ions by red (dark), and -1 
salt ions by gray spheres. Corresponding times 
of snapshots in units of r : (a) 10, (b) 25, (c) 
40, (d) 55, (e) 70, (f) 85. Aggregates of poly- 
electrolytes with scaffold-like structures can be 
seen in the snapshots e and f. Finally, all of 
the polyelectrolytes form a single bundle (last 
configuration is not shown here). 
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Figure 2: Snapshots from evolution of a raft- 
like structure formed by N p = 9 polyelectrolytes 
in the presence of 3-valent salt of concentration 
C3 : i = 0.5. The positive 3-valent salt ions are 
shown by golden (light) spheres and the neg- 
ative salt ions and counterions are not shown 
for clarity. Corresponding times of snapshots in 
units of To: (a) initial configuration before equi- 
libration of the ions, (b) 50, (c) 250, (d) 400, (e) 
550, (f) 700. 
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Figure 3: Potential of mean force between two 
rods as a function of 9 at different values of R 
and salt concentration c^-.i (see the schematic 
configuration of the rods in the inset of the mid- 
dle plot). The legends of the other two plots are 
the same as that of the rightmost plot. 
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Figure 4: Potential of mean force between two 
parallel rods (see the schematic configuration of 
the rods in the inset) function of center- 
to-center separation R for different concentra- 
tions of 3-valent salt. It can be seen that at 
salt concentrations c^i > 0.2 there is attraction 
between two parallel polyelectrolytes. 
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Figure 5: Potential of mean force between two 
perpendicular rods (schematic configuration of 
the rods is shown in the inset) function 
of R at different concentrations of 3-valent salt. 
The size of the circles on the c%.\ = curve 
shows the size of error bars for all our data 
points. Note that the potential is not shown 
for R < 1.6a where the excluded volume inter- 
action takes over and generates a strong repul- 
sion. 
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Figure 6: Preferred angle 6 P and interaction 
type between two charged rigid rods which are 
fixed on two parallel plates of separation R, 
in different regions of Cs-i — R plane (see the 
schematic configuration of the rods in the mid- 
dle plot of Fig. E]). 
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Figure 7: Ensemble average of the interaction 
energy (Lennard- Jones + Coulomb) per particle 
between two rods with their separation fixed at 
R = 2.5a function of 8 with no added salt. 
Inset: potential of mean force between the same 
two rods. The different trends of the two curves 
shows that the entropy of counterions has an 
important role in the orientation of the rods. 
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Figure 8: (a) Evolution time of raft-like struc- 
tures to bundles of parallel polyelectrolytes as a 
function of N p , number of polyelectrolytes form- 
ing these structures. Error bars are obtained 
from averaging over five runs of the simulation 
with different random number sequences. Inset: 
Log-log plot of t e versus N p and a line with the 
slope of 2.7 (dashed line), (b) Activation energy 
extracted from (a) as a function of N p . 
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